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Abstract. We study the response of an ensemble of synchronized phase oscillators to an external harmonic 
perturbation applied to one of the oscillators. Our main goal is to relate the propagation of the perturbation 
signal to the structure of the interaction network underlying the ensemble. The overall response of the 
system is resonant, exhibiting a maximum when the perturbation frequency coincides with the natural 
frequency of the phase oscillators. The individual response, on the other hand, can strongly depend on 
the distance to the place where the perturbation is applied. For small distances on a random network, 
the system behaves as a linear dissipative medium: the perturbation propagates at constant speed, while 
its amplitude decreases exponentially with the distance. For larger distances, the response saturates to 
an almost constant level. These different regimes can be analytically explained in terms of the length 
distribution of the paths that propagate the perturbation signal. We study the extension of these results 
to other interaction patterns, and show that essentially the same phenomena are observed in networks of 
chaotic oscillators. 

PACS. 05.45.Xt Synchronization; coupled oscillators - 05.65.+b Self-organized systems 

1 Introduction prising from mechanical and physico-chemical processes 

|3I4I5| to biological phenomena |6I7I8| . and is reproduced 



Synchronization is a paradigmatic mode of emergent col- 
lective behaviour in ensembles of interacting dynamical el- 
ements |1I2| . It arises in a broad class of real systems, com- 



by a variety of mathematical models. Roughly speaking, 
it consists of some kind of coherent evolution where the 
motions of individual elements are correlated in time. De- 
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pending on the nature of the individual dynamical laws 
and on the interactions, different forms of synchronized 
states are possible. They range from full synchronization, 
where all the elements follow the same orbit in phase 
space, to weakly correlated forms where the ensemble splits 
into almost independent clusters of mutually synchronized 
elements, or where coherence manifests itself in just a few 
state variables or in time averages of suitably chosen quan- 
tities |9I2| . Full synchronization is typical for globally cou- 
pled ensembles of identical elements, where any two ele- 
ments interact with the same strength. When coupling 
is strong enough and represents an attractive interaction, 
the asymptotic state where all elements share the same 
orbit is stable ^Uj. Under certain conditions, stability of 
the fully synchronized state can also be insured for more 
complex interaction patterns, where each pair of elements 
may or may not interact illj . Weaker forms of synchro- 
nization are characteristic of ensembles of non-identical 
dynamical elements. 

A simple but quite useful model for an ensemble of 
coupled dynamical elements is given by a set of N phase 
oscillators, whose individual dynamics in the absence of 
coupling is governed hy (p = oj. The phase <j){t) £ [0, 27r) 
rotates with constant frequency lu. This elementary rep- 
resentation of periodic motion, originally introduced as 
a model for biological oscillations 7 , approximates any 
cyclic dynamics, even in the presence of weak coupling 
13 . As for the interaction pattern, it can be thought of as 
a graph, or network, with one oscillator at each node. The 
graph is characterized by its adjacency matrix = {Jij}. 
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If oscillator i is coupled to oscillator j, i. e. if the phase 
(t) enters the equation of motion of {t) , we have Jij = 
1, and Jij = otherwise. The adjacency matrix is not 
necessarily symmetric and, thus, coupling is not always 
bidirectional. In other words, the interaction network is 
generally a directed graph. The coupled oscillator ensem- 
ble is governed by the equations ^21 

N 

= + fc ^ Jij sin((/)j - 4>i) (1) 

for i = 1, . . . , iV, where k is the coupling strength. The case 
of global coupling, Jij = 1 for all i and j, has been studied 
in the thermodynamical limit by Kuramoto [S] , who found 
that, as the coupling strength grows, the system undergoes 
a transition to a state of frequency synchronization, as first 
predicted by Winfree ^Hl- The transition parameters are 
determined by the distribution of natural frequencies uji. 

For identical oscillators, LUi ^ uj for all i, and for all 
fc > 0, the long-time asymptotic state of a globally cou- 
pled ensemble is full synchronization. More generally, it is 
possible to show that full synchronization is stable when 
the interaction network is regular, i. e. when all oscillators 
are coupled to exactly the same number z of neighbours 
In this situation, = z for all i. 

A fully synchronized oscillator ensemble can be thought 
of as an active medium in a rest-like state. Microscopically, 
this stable state is sustained by the highly coherent collec- 
tive dynamics of the interacting oscillators. A key feature 
characterizing the dynamical properties of the medium is 
determined by its response to an external perturbation. 
How is the synchronized state altered as the perturbation 
signal propagates through the ensemble? Which propaga- 
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tion properties does coupling between oscillators establish 
in the medium? The effect of external forces on ensem- 
bles of interacting dynamical systems has been studied 
in detail for global coupling, both for periodic oscillators 
and chaotic elements |14I15I16I17| . Ordered oscillator ar- 
rays have also been considered The above questions, 
however, are especially significant for more complex in- 
teraction patterns -in particular, for random interaction 
networks- where the non-trivial geometric structure is ex- 
pected to play a relevant role in the propagation process. 
Quite surprisingly, the problem seems to have been ad- 
dressed for the first time only recently |19I2U| . 

In this paper, we present numerical calculations and 
analytical results for the propagation of a perturbation in 
an ensemble of identical phase oscillators, governed by the 
equations 

N 

0i = ^ Ji] sin((/)j - (j>i) + a6ii sm{f2t - (2) 

Without generality loss, the natural frequency of oscilla- 
tors and the coupling strength are fixed to lu = and 
fc = 1, respectively. The external perturbation is repre- 
sented by an additional oscillator of constant frequency 
n, to which a single oscillator in the ensemble -i. e., os- 
cillator 1- is coupled with strength a. To take advantage 
of certain analytical results regarding regular graphs, we 
take a connection network where each oscillator is cou- 
pled to exactly z neighbours. Our main goal is to relate 
the response of the ensemble to the metric properties of 
the interaction network. In particular, we pay attention to 
the dependence on the distance to the node where the per- 
turbation is applied, and on the perturbation frequency. 
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The paper is organized as follows. In the next section, 
we present numerical results for random interaction net- 
works of phase oscillators. We show that the system re- 
sponse is, essentially, a resonance phenomenon, and iden- 
tify two regimes in the dependence on the distance. In 
Sect. 121 we reproduce the numerical results through an 
analytical approach in the limit of small-amplitude pertur- 
bations. We propose an approximation to obtain explicit 
expressions which clarify the role of the network structure 
in the propagation process. In Sect. 0] we extend our re- 
sults to other geometries, including ordered networks with 
uni- and bidirectional interactions and hierarchical struc- 
tures. Moreover, we show that the same propagation prop- 
erties are observed in networks of chaotic oscillators, which 
broadly generalizes our conclusions. Results are summa- 
rized and commented in the final section. 

2 Numerical results 

We have solved Eqs. (0) numerically, for an ensemble of 
~ 10'^ oscillators. We have considered a random regu- 
lar network, with z = 2. The network was constructed by 
chosing at random the z neighbours of each node. Multiple 
directed links between any two nodes were avoided, and 
realizations which produced disconnected networks were 
discarded. Most of the results presented here correspond 
to a perturbation of amplitude a — 10^^ and various fre- 
quencies, typically ranging from Q ^ 10^^ to 10. The in- 
tegration At in our numerical algorithm was chosen such 
that Q-^ > At. 
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The ensemble was prepared in the state of full synchro- 
nization, with (pi — (f)j for all i and j. Before recording the 
evolution of the phases (t>i, an interval much longer than 
f2~^ was left to elapse. After this interval, any transient 
behaviour due to the combined effect of the dissipative 
mechanisms inherent to the coupled oscillator dynamics 
and the external perturbation had relaxed and the system 
had reached a regime of steady evolution. Numerical re- 
sults show that, in this regime, each phase (piit) oscillates 
around the average ^(i) = X^i <^*(^)i seemingly with 
harmonic motion of frequency i7. Our aim was to quanti- 
tatively characterize the departure from full synchroniza- 
tion due to the response to the external perturbation. As 
a measure of this departure for each individual oscillator 
i, we considered the time- averaged mean square deviation 
in (/)-space, defined as 

a^^^[{{q>.-^r)f\ (3) 

where (•) denotes time averages over sufficiently long in- 
tervals. The mean square deviation cr^. is a direct measure 
of the amplitude of the motion of each phase with respect 
to the average (p. In the fully synchronized state, ct^. = 
for all i. 

It turns out that, for a given realization of the interac- 
tion network and a fixed value of the frequency f2, a^. has 
a rather well defined dependence on the distance di from 
oscillator 1, where the external perturbation is applied, to 
oscillator i. The distance di is defined as the number of 
links along the shortest directed path starting at oscilla- 
tor 1 and ending at i. On the other hand, especially for 
large distances di, a^^ may strongly depend on the specific 
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Fig. 1. Mean square deviation from full synchronization as a 
function of the distance from the oscillator where the external 
perturbation is applied, for various values of the frequency fi. 
As a guide to the eye, spline interpolations are shown as curves. 

realization of the interaction network. In view of the well 
defined dependence of the mean square deviation on the 
distance, and for clarity in the notation, from now on we 
drop the index i which identifies individual oscillators. 

Figure n shows results for the mean square deviation 
(70 as a function of the distance d for several frequencies Q 
and a fixed interaction network. In this realization of the 
network, the maximal distance between oscillator 1 and 
any other oscillator is dmax = 15. For each value of d, the 
individual values of have been averaged over all the os- 
cillators at that distance from oscillator 1. It is apparent 
that the mean square deviation exhibits two well differen- 
tiated regimes as a function of the distance. For small d, 
fjfjj decreases exponentially, at a rate that sensibly depends 
on the frequency Q. In this regime, the perturbation is in- 
creasingly damped as it propagated through the system. 
At large distances, on the other hand, is practically in- 
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dependent of d. It is this large-distance value of which, 
typically, shows considerable variations between different 
realizations of the interaction network. The transition be- 
tween the two regimes is mediated by a zone where cr^ at- 
tains a minimum, which is sharper for smaller frequencies. 
Thus, the mean square deviation varies non-monotonically 
with the distance. 

The dependence of cr^ on the frequency fl for fixed 
distance reveals that the response of the oscillator ensem- 
ble to the external perturbation is, essentially, a resonance 
phenomenon. In Fig. [3 we plot as a function of for 
three values of d. For small and large distances, the mean 
square deviation from full synchronization displays a sym- 
metric peak around the natural frequency of the individ- 
ual oscillators {uj = 0). This shows that the response of 
the ensemble is maximal when the external perturbation 
varies with the same frequency of the elementary compo- 
nents of the system. For intermediate distances, however, 
an anomaly appears. While for large frequencies the re- 
sponse decreases as expected, the resonance peak is re- 
placed by a local minimum at = cj. The response is 
now maximal at two symmetric values around the natural 
frequency of the oscillators. This anomaly is directly re- 
lated to the fact that the minimum in for intermediate 
distances (Fig. ^) is sharp for small !? and becomes much 
less distinct as the frequency of the external perturbation 
increases. This effect, combined with the overall decrease 
of tT0 as f2 grows, implies a non-monotonic dependence f2 
at such distances, which results into the appearance of the 
double peak. 
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Fig. 2. Mean square deviation from full synchronization as 
a function of the frequency of the external perturbation, at 
various distances from the node at which the perturbation is 
applied. Curves, added for clarity, are spline interpolations. 

The double- logarithmic plot of Fig.|3reveals the large- 
]7 behaviour of the mean square deviation from full syn- 
chronization. Beyond the resonance zone, decreases 
with the frequency as f2^^ for all distances. Note, how- 
ever, that the value of in the large- regime is the 
same for all d > 0, while for d = -i. e., at oscillator 
1, where the external perturbation is applied- is three 
orders of magnitude larger. 

The fact that the mean square deviation cr^ shows a 
well defined dependence on the distance does not neces- 
sarily imply that the individual motions of oscillators with 
the same value of d are related in any specific way. As a 
matter of fact, being a time average, bears no informa- 
tion about possible correlations between the instantaneous 
state of different oscillators. To detect such correlations we 
have inspected successive snapshots of the set of individ- 
ual phases plotted on the unit circle, i. e. on the plane 
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Fig. 3. The same as in Fig.|21 in log-log scale. Data for d — 
and d = 3 are also shown. Dotted lines have slope —1. 

(cos (p, sin (/)) . Figure 0] shows one of these snapshots for 
the same system of Figs. ^ to 13 It turns out that, actu- 
ally, there is a strong correlation between the positions of 
oscillators with the same value of d. According to their 
distance, they form clusters of gradually decreasing devi- 
ations and growing dispersion. Clusters with d < A are 
so compact that they cannot be resolved into single ele- 
ments in the scale of the main plot of Fig. ^ For larger 
distances, clusters are relatively more disperse, as shown 
in the insert. Clustering in the distribution of phases re- 
veals that, for a given value of d, oscillations around the 
average phase occur coherently. In the next section we 
show analytically that, in fact, these oscillations are very 
approximately in-phase. 

Finally, we have studied the dependence of the re- 
sponse of the system on the perturbation amplitude a. 
Over a wide range, the mean square deviation of individ- 
ual oscillators results to be proportional to the amplitude, 
tj0 (X a. As an illustration. Fig. |S1 shows the ratio a^/a as 
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Fig. 4. Snapshot of individual phases on the plane 
(cos sin 0). Labels indicate the distance to the oscillator 
where the external perturbation is applied. The insert shows 
a close-up of oscillators with d > 3, revealing the fine cluster 
structure for larger distances. 



a function of the distance for a perturbation frequency 
f2 = 0.1. Curves collapse for all a < 2. Only for a > 2 do 
we find significant departures from the small-amplitude 
regime. For a = 3 the response of the whole system is 
relatively increased, especially, for intermediate and large 
distances. Finally, for a = 10 we have an overall saturation 
of the response, and the ratio a^/a decreases. 

In the next section, we show that most of the numer- 
ical results presented here can be analytically explained 
by means of a small-amplitude approximation of Eq. 
Our analytical approach reveals the role of the interaction 
structure in the response of the system to the external 
perturbation. 
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Fig. 5. Normalized mean square deviation, cr^/a, for various 
values of the perturbation amplitude a, as a function of the 
distance to the oscillator where the external perturbation is 
applied. The perturbation frequency is = 0.1. 

3 Analytical results 
3.1 Small-amplitude limit 

The numerical results presented in Fig. suggest that the 
limit of small external perturbation (a ^ 1) bears signif- 
icant information on a wide range of values for the am- 
plitude a. It is therefore worthwhile to study this limit 
analytically, representing the instantaneous state of the 
oscillator ensemble as a perturbation of order a to full 
synchronization. Introducing, as in the previous section, 
the average phase (t>{t) — N^^ J2i 4'i{t): we write the phase 
of oscillator i as a perturbation of order a to 0, 

(i3,{t)^4,{t)+ajp,{t), (4) 

with = 0. In turn, the average is expected to vary 

around a constant phase ^o, with fluctuations of ampli- 
tude proportional to a: 

^{t) = dpo + a^{t). (5) 
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Without generality loss, we take = 0. 

Replacing Eqs. I@J and (O in (01, and expanding to 
the first order in the perturbation amplitude a, we get 

^ = 1 ^ ~^^) + ]^ exp(ir2t) (6) 

for the average phase deviation, and 

= -<i> + ^ Jij {ipj - i/'O + Sii ex.p(int) (7) 

3 

for the individual deviations. For simplicity in the mathe- 
matical treatment, we have replaced shi(f2t) by exp(il?i). 
Focusing on the case where each oscillator is coupled to 
exactly z neighbours, the above equations can be simpli- 
fied taking into account that Jiji'^Pj ~ ''Pi) — —zipi + 
and 'Zij Jijii^j -^i) = Jiji^j- 
Note that the time derivative of the average phase de- 
viation <P enters the equation of motion for the individual 
deviations ipi, Eq. lO, as a kind of external force act- 
ing homogeneously over the whole ensemble. According 
to Eq. ©, this effective force is of order N^^. As we show 
later, it dominates the response of oscillators at large dis- 
tances, where the effect of the perturbation signal propa- 
gated through the network is lower. 

Equation ||7J| admits solutions of the form ipi{t) = 
exp(ii7t), corresponding to steady harmonic motion at 
the frequency of the external perturbation. These steady 
solutions are expected to represent the motion once tran- 
sients have elapsed. The complex amplitudes Ai satisfy a 
set of linear equations, which can be cast in matrix form 
as 

CA = b, (8) 
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with A — (Ai, ^2, • ■ ■ , J^n)- The elements of the matrix 
£ and of the vector b are, respectively, 

k 

and 

h = ~ ^- (10) 

For a given interaction network, the amplitudes A = £^^b 
can be found, for instance, numerically. Note that the 
amplitude modulus \Ai\ is a direct measure of the mean 
square deviation cr^. , introduced in Sect.|21to quantita- 
tively characterize the departure from full synchroniza- 
tion in the mmierical realization of our system. In fact, 
for harmonic motion, a^^ — a|^j|/V2. The phase ipi of 
the complex amplitude Ai = exp(i(pi) measures the 
phase shift of the oscillations of i)i[t) with respect to the 
external perturbation. 

Figure El shows the amplitude moduli \Ai\ obtained 
from the numerical solution of Eq. © , for a ensemble of 
TV = 10'^ oscillators with the same interaction network as 
in the results presented in Sect. [3 The values of \Ai\ are 
averaged over all the oscillators at a given distance d from 
oscillator 1, and the subindex i is accordingly dropped. 
Results are presented for several values of the frequency 
fl. Dotted lines are spline approximations of the numerical 
data for already presented in Fig. ^ multiplied by a 
factor \/2/a. The agreement with the solution of Eq. ||HJ) 
is excellent. 

The phase shifts ipi corresponding to the same solu- 
tions of Eq. (jSJ are shown in Fig. [T] as a function of the 
distance and for several values of the frequency fl. Data for 
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Fig. 6. Amplitude moduli from the numerical solution of 
Eq. 0, for the same system as in Fig.Q Dotted lines are spline 
approximations for \/2a^/a, for the values of the mean square 
deviation presented in Fig. Full straight lines show the 
analytical prediction for the slope at small distances, Eq. 11611 . 

n — 1 are also shown. As it may be expected, phase shifts 
are always negative, indicating a delay in the response of 
the system to the external perturbation. Coinciding with 
the regime where the amplitude moduli \Ai\ decrease ex- 
ponentially, we find a zone where phase shifts vary linearly 
with distance. Namely, the phase shift between oscillators 
whose distances to oscillator 1 differ by one is a constant 
Aip. This unitary phase shift, which gives the slope of 
as a function of the distance, depends on the frequency 
]7. In this zone, individual phases vary with distance as 
V'i oc exp[i((iiZi(^ + Qt)]. Thus, the perturbation propa- 
gates through the system at constant speed |i7/zii^|. For 
larger distances, this linear regime breaks down, and the 
variation of with d tends to be much less pronounced. In 
the transition between both regimes, however, the phase 
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Fig. 7. Phase shift as a function of distance, from the numer- 
ical solution of Eq. Q and for the same system as in Fig. Q 
Straight lines show the analytical prediction for the slope at 
small distances, Eq. (1161 . 

shift varies rather abruptly for smaU frequencies, while it 
develops a minimum for larger Q. 

The results presented in Figs. El and [3 where obtained 
from the numerical solution of the small-amplitude limit 
equation 0. A more explicit solution can be obtained 
from a suitable approximation of Eq. , taking into ac- 
count specific mathematical properties of the matrix £, as 
we show in the following section. 

3.2 Approximate solution 

According to Eq. @, the matrix C can be written as 

C^{z + iQ)I-J, (11) 

where T is the N x N identity matrix, and the elements 
of J are 



with — Jkj ■ Our approximation to the solution of 
Eq. ^ is based on the following remarks. 

(i) Due to the fact that J^j Jij — ^ for h the eigen- 
values of the adjacency matrix are all less than or equal 
to z in modulus Moreover, the eigenvalues of J are 
the same as those of J^. In fact, if v = (vi,V2, ■ ■ ■ ,vn) 
is an eigenvector of J^, then v = {vi,V2, ■ ■ ■ ,vn), with 
V'i — Vi — iV^^ J2k'^k, is an eigenvector of J' with the 
same eigenvalue. This implies that, for f2 ^ 0, the inverse 
of the matrix C can be expanded as 



(13) 



(z + in) ^ (z + i/?)™+ 

^ ^ m— 

because all the eigenvalues of the matrix (z + \fl)^^J are 
less than unity in modulus. 

(ii) While = z for all i, J^j = varies with 

j. Note that is the number of links starting at j, and 
thus gives the number of oscillators which are coupled to 
oscillator j. For any realization of the interaction network, 
however, the average value of over the whole ensemble 
is always the same, J^j — ^jk ^kj — z. This 
suggests that, as an approximation to Eq. (|12|l avoiding 
the explicit calculation of ^j, WG 00/11 tj0,kG Jij "~ Jij ' 
z/N. More generally, it is possible to show that this same 
approximation yields, for the powers of J^, 



T(m) _ j(m) _ 



(14) 



Jij — Jij jy^J' 



(12) 



where jiT^-' and Jl^-' are elements of the matrices J"^^ and 
''J '■J 

>7™, respectively. 

Combining Eqs. (|13|l and H14|) . the approximate form 
of matrix C^^ can be applied to the vector b in the right- 
hand side of Eq. © to give the following approximation 
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for the amplitudes: length between any two oscillators is, at most, of order 



m— 



j{m) _ _Z_ 

^1 N 



z/N (which we assume to be a small parameter, as in our 
numerical analysis). As a result, for most oscillators at a 

oo j(m) . 

E7 TT^T — PT + (15) Small distance from oscillator 1, we have J^f^'^ — 1. More- 



oo 

' ^ (z + if2)"+i 



^ (z + ii7)"+i Nfl 

m— ^ ' 

In this approximation, the effect of the average devia- 
tion from full synchronization -which, as we discussed in 
Sect. 13. ll can be interpreted as an external effective force 
acting over all the oscillators with the same intensity- is 
represented by the terms of order A^^^. 

It is interesting that the approximation H15|) involves 
explicitly the matrix elements of the powers of J . The ma- 
trices J™ {m — 1,2,...) bear information about the met- 
ric structure of the interaction network. Specifically, the 
element J^j™'' equals the total number of directed paths of 
length m starting at node j and ending at node i |21| . In 
other words, ■/f™' gives the number of different ways of 
reaching node i from node j in exactly m steps along di- 
rected links. Equation IjlSfl shows that the response of any 
individual oscillator to the external perturbation, mea- Within this approximation, the small-distance exponen 



over, the total number of nodes with small di, n{di) ^ z'^' , 
is also small as compared with the system size N. This 
implies that the possibility that an oscillator at a small 
distance di is also connected by a path of length slightly 
larger than di can be neglected. Consequently, for oscilla- 
tors at small distances from the node at which the pertur- 
bation is applied, the sum in the second line of Eq. (|15|l 
is dominated by the term with m = di. This dominance 
is enhanced for large \z + iJ7|, because successive terms 
in the sum are weighted by increasing inverse powers of 
that number. If the system is large enough we can drop 
the last term in the second line of Eq. IjlSI) . and write 
A, w {z + in)-'^^-\ i. e. 



A, w {z'' + n^) exp 



-i{di + l)tan ^ — 
z 



(16) 



sured by the amplitude Ai, is directly related to the num- tial dependence of the amplitude modulus, |^,;| « (z^ + 

ber of paths though which the perturbation signal can i?^) 2— is apparent. Straight lines in Fig. ^ show the 

flow from oscillator 1. Note that jf™' = for m < di excellent agreement between the predicted slope of \Ai\ 

and Jj^™'' = 1 for m = di. For oscillator i, therefore, the and our numerical results. Equation (|16|l also explains 

first contribution to the sum in the second line of Eq. IjlSfl the linear dependence of the phase shift ipi with the dis- 

comes from the term with m = di. For m > di, J^^™"* tance. Specifically, it predicts a unitary phase shift Aip — 

different from zero if at less one path of length m starts — tan^^(r2/z). Straight lines in Fig. [7|stand for this pre- 

at oscillator 1 and ends at i. diction. 

For nodes at small distances from oscillator 1 , there is Two effects contribute to break down the small-dist- 

typically only one path of length di from 1 to i. In fact, ance approximation (|16|l . First, as discussed above, we 

the probability of having more than one path of short expect that this approximation does not hold beyond dis- 
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tances where z'^' ^ N, i. e. di ^ log N/ log z. At larger This large-distance approximation for the complex ampli- 

distances, in fact, it is not true that only one path con- tude Ai is more clearly analyzed for limit values of the 

tributes to the propagation of the perturbation signal from frequency f2. For small J7, specifically for dif2/z <ti 1, the 

the node at which it is applied. The second effect has to amplitude modulus is 



do with the relative magnitude of the first non-zero term ^ j d^Q'^ 
in the sum of the last line of Eq. H15|) and the term of 

order . For sufficiently large distances, the latter can- As found in our numerical results, Figs. (HJ and ©, the 

not be neglected with respect to the former. If < z, the response grows with d, for large distances. Combined with 

two terms become comparable for z'^'+i - NQ, while if the decrease observed for small distances, this explains the 

z < I? they are similar for fi'^- - N. It turns out that, existence of an intermediate minimum in both \Ai\ andcr^. 

in both limits, this second effect acts at distances smaller The amplitude phase 

than the first one. This does not imply, however, that the i z / 1 \ 
first effect plays no role in determining the response of the 

system at large distances. exhiihiis a more complicated functional dependence with 

It can be shown that, for osciUators at large distances distance. In the opposite Umit of large frequencies, the 

from oscillator 1 , the matrix clement in the sum of Eq. m l^'-g^-distance approximation ^ is dominated by the last 

can be accurately approximated as J^r^ = Jo^" for all t^™' ^» ~ We recall that this term stands for the 

m > d„ while - for m < d,. In fact, we can argue contribution of the average deviation from full synchro- 

that the number of paths of length m ending at a given nization to the individual motion of oscillators. The am- 

oscillator i scales as z™ for large m, by noticing that this P^it^^^^ modulus becomes independent of the distance, as 

number is z times the number of paths of length m - 1 ^^^^^ in the results of Fig. E Its phase is always 7r/2 -or, 

ending at the oscillators to which i is coupled. The precise '^o^''^ generally, ^/2 + 2fc7r, with k an integer. In the results 

value of the prefactor Jo depends on N and on the specific ^^^^ l^^^S^ ^ ^^^"^^^ Fig.Q we have ^ « 7r/2-2^ « -4.71. 

realization of the network, but is independent of the dis- ^o^^i^^S ^ack to the full form of our approximation 

tance d,. For the network corresponding to the numerical f""" ^^"^ amplitude, Eq. let us finally point out that, 

results presented above, for instance, Jo « 9.91 x 10"^. f""" sufficiently large frequency f2 and irrespectively of the 

Replacing the Ansatz for J^r^ in Eq. m and per- distance d„ the dominant terms arc of order f2-\ To this 

forming the summation, we get O'^^e'^' contributes its first term, m = 0. Since 



-Jo [ f H 2" ) cxp [ —idi tan 



= X, for J7 ^ oo we get 

z ) ' N ' 

(17) "^^T^r^^^^iV 



A.^^(-6a + l]. (20) 
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This result explains the decay as Sl^^ in the tails of the 
resonance peaks, displayed in Fig.O It also shows that the 
large-frequency response of oscillator 1 is times larger 
than that of any other oscillator, as illustrated in the same 
figure. Moreover, we find that the phase shift of oscillator 
1 is ifi = — 7r/2 while any other oscillator is dephased by 
7r/2 with respect to the external perturbation. 

4 Extensions 

4.1 Other regular network structures 

It is important to remark that the analytical approach 
presented in Sect. |31 so far applied to random networks, 
is valid for a large class of interaction patterns. In fact, 
the only condition imposed on the adjacency matrix 
to obtain the results of Sect. 13.11 is that J^j Jij — ^ foi' 
all i, while the approximation of Sect. 13.21 requires that 

Jjj « z for all j. These conditions imply, respectively, 
that each oscillator is coupled to exactly z neighbours and 
that, in turn, the number of oscillators coupled to each 
oscillator is also approximately constant. Under such con- 
ditions, our approach can be used to evaluate the response 
of an ensemble with any interaction pattern. In this sec- 
tion, we illustrate this fact with a few cases that admit to 
be worked out explicitly. 

Consider first a linear array of N oscillators with pe- 
riodic boundary conditions, where each oscillator is cou- 
pled to its nearest neighbour to the left (z = 1). We as- 
sume that oscillators are numbered from left to right in 
the natural order. For this directed ring, we have 1 if 
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i = {j + 1) mod N, and Jy = otherwise. Consequently, 
Sfc '^kj — 1 for all j, and the approximation of Sect. 13.21 
is exact. The relevant elements of J^™ are j/™'' = 1 if 
m = di + kN {k = 0, 1, 2, . . .), and J^^-' = otherwise. In 
fact, the only paths that join oscillator 1 with oscillator i, 
at distance di — i — 1, are those of length di plus an integer 
number k of turns around the ring. The calculation of the 
amplitudes yields 

(1 + in)-'^'-'^ i 

^ " i-{i + in)-N + Jw- ^^^^ 

For n <^ 1, but with Nf] ^ 1, the amplitude at essen- 
tially all distances {di < N) is dominated by the first 
term of Eg. 1211 In this regime, \Ai\ decreases exponen- 
tially with di and the phase shift ipi varies linearly, with 
Aip = — tan^^ i7. In the opposite limit of large frequen- 
cies, ^ 1, the amplitude decays as ]7~^. For oscillator 
1 {di = 0), we have Ai w — i/i7, while for any other os- 
cillator {di > 0), Ai « i/NQ. We stress the qualitative 
similarity between these results and those obtained for 
random networks. 

Note that in the limit N —i- oo, the perturbation can- 
not attain the oscillators to the left of oscillator 1 , an there 
is only one path (of length m = di) for the signal to reach 
the oscillators to its right. In this limit, Ai = (l-l-ii?)^'''^^ 
at all distances, and the regime of exponential decay and 
linear phase shift extends over the whole system to the 
right of the node where the perturbation is applied. 

Consider now a linear array -which, for simplicity, we 
treat in the limit N oo- where each oscillator is cou- 
pled to its two first neighbours {z — 2). Again, the ap- 
proximation of Sect. 13.21 is exact. With this bidirectional 



Damian H. Zanette: Disturbing synchronization; Propagation of perturbations in networks of coupled oscillators 13 

coupling, there are infinitely many directed paths joining contributes to the propagation of the perturbation signal 

any two nodes in the network. A path of length m be- towards each oscillator. For the bidirectional linear array, 

tween oscillator 1 and oscillator i consists of to-|_ steps to in contrast, the total number of paths ■^ii"'' grows 

the right and to_ steps to the left, with m+ + m_ ~ m exponentially with di. This growth compensates partially 

and |m+ — m_ \ = di. The number of such paths is the exponentially decreasing contributions from successive 







( m \ 















=1 1 = 1 1 (22) 



path lengths [cf. Eq. ifT^ ]. 

In the limit of large frequency, on the other hand, we 



IS essen- 



■ r J J ^ ^ J ^ u ^x, AA ^ n ^-1 find w (2 + iny^'-^, i. e. \AA - Q''^'-^. It is 

if m and di [m > di) are both even or odd, and other- ' v / ' i 'i 

m r ,1 i- • TTTT -i • tial to tliis result, however, that -assuming an infinitely 

wise, io perform the summation m Eq. i|15|i . it is conve- ' ' o j 

. ., 7,r, J ni rni- large system- we have neglected the term of order A^~^ in 

nient to write m = di + Zq and sum over (7 = 0,1,.... i liis o j o 

yj^y^ Eq. H15|l . For any finite size, if N is fixed, the limit of large 



-1 




frequencies is dominated by this term, and \Ai\ ~ ]7 



Let us finally consider a interaction network with a hi- 



— ' (2 + iJ7)'^i+29+i erarchical structure, in the form of a directed tree where 

q=0 ^ ' 



T? \di±\ di 



4 



(2+il7) 



nodes are distributed into successive layers. An oscilla- 



(23) 



(2 + i]7)'''+i g^-j- ^ given layer is coupled to only one oscillator in the 

where 2Fi{a,b;c;z) is the hypergeometric function 24. layer immediately above, so that z = 1, and the uppermost 

This exact result can be approximately analyzed for lim- layer consists of a single oscillator, where the external per- 

iting values of f2. For i7 ^ 1, the amplitude modulus turbation is applied. Thus, the perturbation propagates 

decreases exponentially with the distance, as \Ai\ oc (1 + downwards through the hierarchy. If layers are labeled in 

y/2n)-'^-/'^, while the phase shift varies linearly with di, the natural order starting by / = at the uppermost layer. 



with slope Aip ~ — tan^^(l + . Note that the the distance of any oscillator to the node where the per- 

exponential decrease of \Ai\ is much slower than in the turbation is applied coincides with the label of its layer, 

case of random networks with the same z. For random Namely, oscillators at level / = 1 have di = 1, at level 

networks, in fact, we have found that -in the limit of 1 = 2 they have d,; = 2, and so on. There is only one path, 

small frequency and in the regime of exponential decay- of length di, through which the perturbation can reach os- 

\Ai \ ~ z"''*"^. Here, on the other hand, the decay of \Ai\ dilator i. Therefore, only one term contributes to the sum 

becomes increasingly slower as J7 ^ 0. This important in Eq. H15|) . Now, it is important to note that the approx- 

difference is a consequence of the fact that, in a random imation (ii) introduced in Sect. 13.21 is no longer suitable, 

network and at small distances, essentially only one path While each oscillator is coupled to only one neighbour. 
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the number of oscillators coupled to oscillator J, X^fc '^fcji 
is typically larger than one. In a regular hierarchical struc- 
ture, in fact, we fix J^j — z' > 1 for all j . In this sit- 
uation, Eq. p5(l changes in such a way that the last term 
becomes multiplied by a factor z' /z. The amplitude is 

A. ^(1 + 1/2)--- + ^^. (24) 

Not unexpectedly, this result is similar to that for a uni- 
directional linear array, Eq. H21|) . In the first term, which 
stands for the exponential-decay regime, the main differ- 
ence corresponds to a factor which, in the case of the array, 
takes into account that the perturbation signal reaches an 
oscillator at each turn around the array. In the second 
term, the two results differ precisely in the factor z' / z. 
This effect enhances the response of oscillators with di > 
at large frequencies, where \Ai\ w z' /zNQ. 

4.2 Non-regular random networks 

So far, our numerical and theoretical analyses have dealt 
with regular networks, where all oscillators are coupled to 
exactly the same number z of neighbours. While the an- 
alytical approach cannot be extended to the case of more 
general structures, it is worthwhile to show that our results 
hold -at least, qualitatively- for non-regular random net- 
works. For sufhciently large networks, with a well defined 
average number of neighbours per site, it is in fact ex- 
pected that statistical quantities such as the mean square 
deviation cr^. are essentially not sensible to the regularity 
of the interaction pattern. 
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We consider non-regular random networks of two types. 
In the first type (random I) the number of neighbours Zi 
of each site i is chosen to be 1, 2 or 3 with equal prob- 
ability 1/3. The average number of neighbours is thus 
z = 2, which makes it possible to compare with our re- 
sults for regular networks with z = 2. Once Zi has been 
defined, the neighbours of site i are chosen at random from 
the whole system, avoiding self-connections and multiple 
connections. In the second type of non-regular random 
networks (random II), the number of neighbours of each 
site is drawn from the discrete probability distribution 
p{z) = 2^^ (z = 1,2, . . .), which also insures z = 2 but 
with a much wider dispersion. Our numerical calculations 
were run for a system of = 10^ oscillators, with a per- 
turbation amplitude a — 10^"^. 

Figure |H1 shows results for the mean square deviation 
from full synchronization a^, with the two types of non- 
regular random networks. For the sake of comparison with 
the case of regular networks, the perturbation frequencies 
fl are the same as in Fig. ^ Moreover, specific realiza- 
tions of the networks with the same maximal distance to 
the perturbed oscillator, (i,„ax — 15, were selected. The 
vertical axes cover also the same range. We verify at once 
that the main features in the dependence of the mean 
square deviation on the distance found for regular inter- 
action patterns are also present in non- regular networks. 
As it may have been expected, quantitative differences are 
more important for small frequencies, i.e. near the reso- 
nance. There, the oscillator network is more sensible to the 
perturbation and, arguably, its detailed structure plays a 
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Fig. 8. Mean square deviation from full synchronization in 
non-regular random oscillator networks of the two types de- 
scribed in the text (random I and II) , for various values of the 
perturbation frequency Q. Symbols are as in Fig. 

more noticeable role in determining its response. For larger 
frequencies, the values of cr^ become increasingly indistin- 
guishable from those obtained for regular networks. 

4.3 Chaotic oscillators 

An important question regarding the generality of the re- 
sults presented so far is whether they apply to ensembles 
of coupled elements whose individual dynamics are not 
simple phase oscillations. While we may argue that any 
cyclic behaviour, even in the presence of external forces. 
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can be approximately described by a periodic phase os- 
cillator [Sj , the question remains open for chaotic coupled 
dynamical systems. To address this problem we have con- 
sidered an ensemble of Rossler oscillators, described by 
the equations 

ii = -Vi - Zi + k J2j Jij [xj - Xi) + aSii sin Qt 
Vi ^ Xi+ 0.2yi + k J2j Ji] iVj - Ut) 

Zi 0.2 + Zi{xi - c) + k Y,j Jij{zj - Zi). 
The parameter c controls the nature of the oscillations; 
for c = 4.46 they are chaotic |221- As in Sects. |21and|21 
we choose the adjacency matrix such that = z for 

all i. Each Rossler oscillator is thus coupled to exactly z 
neighbours. Chaotic systems can be fully synchronized if 
the coupling intensity k is larger than a certain threshold 
value, related to the Lyapunov exponent of the individual 
dynamics |2I11| . For the above value of c, a coupling in- 
tensity k — 0.2 insures that full synchronization is stable. 
In our ensemble of Rossler oscillators, the external pertur- 
bation acts on just one of the coordinates of oscillator 1, 
namely, on xi{t). 

The numerical results presented below have been ob- 
tained for an ensemble of iV = IQ-^ Rossler oscillators, 
with z = 2 and the parameters quoted in the preceding 
paragraph. For the sake of comparison, the interaction 
network is the same as in our study of phase oscillators 
(Sects. 121 and O. The amplitude of the external pertur- 
bation is a = 10^'^. As a characterization of the response 
of the system, we have used a natural extension of the 
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mean square deviation from full synchronization defined 
in Eq. JJ)) for phase oscillators, given by 

ar. = ((|r,-fp))^/^ (26) 

with Ti = (xi,yi,Zi) and f ~ -^^^X^i'**- Time averages, 
indicated as (•), are performed over sufficiently long in- 
tervals, after transients have been left to elapse. Figure |51 
shows numerical results for the mean square deviation as a 
function of the distance, and for various perturbation fre- 
quencies. Individual values of a^^ for Rossler oscillators at 
a given distance from oscillator 1 have been averaged, and 
the index i has been dropped accordingly. We see that the 
overall behaviour of as a function of d is similar to that 
found for ensembles of phase oscillators (cf. Fig. with 
rapid decrease for small distances and smooth growth for 
large distances. In all cases, the two regimes are separated 
by a well defined minimum. The exponential character of 
the decrease for small d and the saturation of for large 
d are, however, much less clear than for phase oscillators. 

The results of Fig.^clearly show that the mean square 
deviation varies non-monotonically with the perturbation 
frequency. The dependence of with f2 is expected to 
reveal the resonance nature of the system response. The 
natural frequency lu of individual Rossler oscillators can 
be defined, in the chaotic regime, in terms of the average 
period T of the chaotic oscillations. Each period is deter- 
mined, for instance, as the time needed for an oscillator to 
cross the plane y = in the subspace a; < 0. These times 
are then averaged over a large number of oscillations, and 
the frequency is calculated a.s lu — 2tt/T . For c = 4.46, we 
find uj w 1.08. 




Fig. 9. Mean square deviation from full synchronization as a 
function of the distance for an ensemble of Rossler oscillators, 
for various values of the frequency Q. Spline interpolations are 
shown as curves. 



In Fig. 1101 we present numerical results for the mean 
square deviation from full synchronization as a function of 
J7, for several distances. For d = 0, we find the expected 
resonance maximum at J7 w Interestingly enough, there 
is an additional local maximum at J? w 2uj, correspond- 
ing to a harmonic resonance induced by non- linear effects. 
The presence of this extra peak is consistent with the fact 
that higher-harmonic components are very relevant contri- 
butions to the chaotic motion of individual Rossler oscil- 
lators pS; . For larger distances, the resonance maximum 
is replaced by a double peak, as we have found for phase 
oscillators (Fig. with a relative mininmm at i? w w 
and two lateral maxima. At the site of the harmonic res- 
onance we find the same structure. In contrast with the 
case of phase oscillators, however, the double peak persists 
at large distances, with better defined maxima as d grows. 
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Fig. 10. Mean square deviation from full synchronization as a 
function of perturbation frequency for an ensemble of Rossler 
oscillators, at various distances from the node where the per- 
turbation is applied. Spline interpolations are shown as curves. 

We have also verified that, as for phase oscillators, the 
individual motions of Rossler oscillators at the same dis- 
tance from oscillator 1 are in-phase. These coherent dy- 
namics give rise to a clustered distribution in r-space, and 
a snapshot of the ensemble in that space produces a pic- 
ture qualitatively very similar to Fig. 01 

5 Discussion and conclusion 

In this paper, we have studied the response of an ensem- 
ble of fully synchronized oscillators to an external pertur- 
bation. The perturbation is represented as an additional 
oscillator, evolving autonomously with a fixed frequency. 
One of the oscillators of the ensemble is coupled to this 
additional element. The perturbation propagates through 
the system due to the coupling between oscillators. In the 
absence of the external action, this interaction sustains 
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the state of full synchronization. The system can thus be 
thought of as an active extended medium with a highly 
coherent rest state -full synchronization- whose response 
to the external perturbation is driven by the collective 
dynamics of the interacting oscillators. Such response, in 
fact, provides a characterization of the collective dynam- 
ics. 

Our study was mainly focused on ensembles of identi- 
cal phase oscillators. In its usual formulation, Kuramoto's 
model considers global coupling, where interactions are 
identical for all oscillator pairs. In this situation, all oscil- 
lators are mutually equivalent, and the response to the ex- 
ternal perturbation -other than on the oscillator where the 
perturbation is applied- is homogeneous over the whole 
system. Therefore, we have considered more complex in- 
teraction patterns, introducing interaction networks which 
allow for a non-trivial distribution of distances between 
oscillators. Interactions were not necessarily bidirectional, 
so that coupling was not always symmetric. To take ad- 
vantage of certain analytical results on the stability of 
full synchronization we have considered regular net- 
works, were all oscillators are coupled to the same number 
of neighbours. Numerical results obtained show, however, 
that this choice does not represent a strong restriction on 
the interaction network. The number of neighbours of each 
oscillator and the frequency of the external perturbation 
are the main parameters that control the response of the 
system. 

For moderate values of its amplitude, the external per- 
turbation induces oscillations around the state of full syn- 
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chronization. Our main conclusion, first found by numer- 
ical means, is that the response of each individual oscilla- 
tor exhibits a clear dependence on the distance from the 
node where the perturbation is applied. For random inter- 
action networks, this dependence shows two well defined 
regimes. At small distances, the amplitude of the individ- 
ual oscillations decreases exponentially with the distance. 
Meanwhile, the phase shift of these oscillations with re- 
spect to the external perturbation, which measures the 
delay of the individual response, varies linearly with the 
distance. In this regime, thus, the perturbation propagates 
through the system at constant velocity, and is progres- 
sively damped at a rate proportional to its own amplitude. 
In other words, the system behaves as a linear dissipa- 
tive medium. For large distances, on the other hand, the 
individual response saturates and the dependence of the 
amplitude and the phase shift with the distance becomes 
much smoother. The relative extension of the two regimes, 
as well as the rates of attenuation and dephasing of the 
perturbation signal in the linear regime, depend on the 
frequency of the external perturbation and on the num- 
ber of neighbours of each node. 

The fact that the phase shift of individual oscillations 
with respect to the external perturbation is defined by the 
distance to the node where the perturbation is applied, 
implies that all the oscillators at a given distance respond 
to the perturbation coherently. Since the amplitudes of 
their oscillations are also similar, the "spatial" distribu- 
tion of the ensemble -i.e. the distribution in the relevant 
one-particle state space- becomes clustered. Especially in 
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the small-distance regime, where the individual response 
strongly depends on the distance, all the oscillators at a 
given distance form a compact cluster with coherent oscil- 
latory motion. Thus, the external perturbation induces a 
"spatial" organization associated with the internal struc- 
ture of the interaction network. 

The overall response of the system to the external per- 
turbation is maximal when the perturbation frequency is 
equal to the natural frequency of the oscillators. This not 
unexpected resonance phenomenon is revealed by the pres- 
ence of a peak in the amplitude of individual oscillations 
as a function of the perturbation frequency. Far from the 
peak, the amplitude decreases as the inverse of this fre- 
quency. For oscillators at intermediate distances, however, 
an anomaly in the response appears. The resonance peak 
is replaced by a double peak, with a minimum at the res- 
onance frequency and two lateral maxima. As we discuss 
below, this effect can be interpreted as an interference 
phenomenon. 

Our numerical results are well reproduced by an ana- 
lytical approach based on a linear approximation for small 
perturbation amplitudes. This approach is able to dis- 
cern between the roles of diflFerent contributions to the 
response of the system. Two complementary aspects are 
worth mentioning. First, we have found that the existence 
of two regimes in the distance dependence is directly as- 
sociated with the distribution of the number of paths in 
the interaction network. The perturbation signal reaches 
oscillators at small distances essentially through only one 
path. The dissipative mechanisms inherent to the dynam- 
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ics of the oscillator ensemble progressively attenuates the 
signal, which decays exponentially with the distance. For 
oscillators at large distances, on the other hand, the num- 
ber of available paths grows exponentially, at a rate that 
-at least, for small perturbation frequencies- is similar to 
the rate of exponential decay of the signal. These two par- 
tially compensating effects determine that the variation of 
the response with the distance is much smoother for large 
distances. 

The second aspect has to do with the fact that the 
external perturbation affects individual oscillators in two 
ways. Besides the propagation of the signal through the 
network, which acts on the oscillators with different inten- 
sities depending on their distance to the node where the 
perturbation is applied, there is an overall contribution 
originated by the average motion of the whole ensemble, 
which affects all oscillators with the same intensity. This 
global contribution is inversely proportional to the system 
size, and therefore could be generally neglected for suffi- 
ciently large systems. However, it does play an important 
role in determining the system response at distances where 
the propagating signal has been strongly damped. It also 
determines the response at high perturbation frequencies, 
where the propagation mechanism is very ineffective, es- 
pecially, at large distances. Our analysis shows that the 
overall contribution of the average motion of the ensemble 
and the local contribution of the propagated signal have 
opposite signs. In other words, their oscillation phases dif- 
fer by TT, in such a way that, if their amplitudes are similar, 
a phenomenon of destructive interference takes place. In 
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our random networks this happens, precisely, at the tran- 
sition between the two regimes discussed above, when the 
amplitude response has decreased by a factor of about the 
inverse of the system size. The minimum in the amplitude 
at those intermediate distances can thus be interpreted as 
the result of the destructive interference between the two 
contributions that affect individual motions. 

Generally, we may expect that interference phenom- 
ena play an important role in the dynamics of oscillator 
networks. This is due to the fact that, as we have seen, 
the oscillatory signal changes its phase as it propagates 
through the network. The signal can reach a given oscilla- 
tor through different paths, and thus with different phases. 
The sum of all those contributions will depend not only 
on their amplitudes but also on their relative dephasing, 
likely giving rise to interference. Prom this perspective, the 
saturation of the response at large distances, were contri- 
butions from many different paths are acting, could be 
interpreted as a phenomenon of constructive interference 
that breaks down the regime of exponential decay. 

The applicability of our analytical approach is not re- 
stricted to random networks. We have shown how our 
main results extend to regular and hierarchical arrays of 
phase oscillators. Numerical calculations show, moreover, 
that the same results are qualitatively reproduced in non- 
regular random interaction patterns. The most important 
extension considered here, however, has consisted of re- 
placing phase oscillators by chaotic Rossler oscillators. 
The response of ensembles of these chaotic elements to the 
external perturbation was analyzed by numerical means. 
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